The present invention relates to a method and an apparatus for image reconstruction in nuclear medicine imaging and, more particularly, but not exclusively to image reconstruction in nuclear medicine imaging using gating techniques.
Radionuclide imaging aims at obtaining an image of a radioactively labeled substance, that is, a radiopharmaceutical, within the body, following administration, generally, by injection. The substance is chosen so as to be picked up by active pathologies to a different extent from the amount picked up by the surrounding, healthy tissue in consequence; the pathologies are operative as radioactive-emission sources and may be detected by radioactive-emission imaging. Pathology may appear as a concentrated source of high radiation, that is, a hot region, as may be associated with a tumor, or as a region of low-level radiation, which is nonetheless above the background level, as may be associated with carcinoma.
A reversed situation is similarly possible. Dead tissue has practically no pick up of radiopharmaceuticals, and is thus operative as a cold region.
The mechanism of localization of a radiopharmaceutical in a particular organ of interest depends on various processes in the organ of interest such as antigen-antibody reactions, physical trapping of particles, receptor site binding, removal of intentionally damaged cells from circulation, and transport of a chemical species across a cell membrane and into the cell by a normally operative metabolic process. A summary of the mechanisms of localization by radiopharmaceuticals is found in http://www.lunis.luc.edu/nucmed/tutorial/radpharm/i.htm.
The particular choice of a radionuclide for labeling antibodies depends upon the chemistry of the labeling procedure and the isotope nuclear properties, such as the number of gamma rays emitted, their respective energies, the emission of other particles such as beta or positrons, the isotope half-life, and the decay scheme.
In PET imaging, positron-emitting radioisotopes are used for labeling, and the imaging camera detects coincidence photons, the gamma pair of 0.511 Mev, traveling in opposite directions. Each one of the coincident detections defines a line of sight, along which annihilation takes place. As such, PET imaging collects emission events, which occurred in an imaginary tubular section enclosed by the PET detectors. A gold standard for PET imaging is PET NH3 rest myocardial perfusion imaging with N-13-ammonia (NH3), at a dose level of 740 MBq, with attenuation correction. Yet, since the annihilation gamma is of 0.511 Mev, regardless of the radioisotope, PET imaging does not provide spectral information, and does not differentiate between radioisotopes.
In SPECT imaging, primarily gamma emitting radioisotopes are used for labeling, and the imaging camera is designed to detect the actual gamma emission, generally, in an energy range of approximately 11-511 KeV. Generally, each detecting unit, which represents a single image pixel, has a collimator that defines the solid angle from which radioactive emission events may be detected.
Because PET imaging collects emission events, in the imaginary tubular section enclosed by the PET detectors, while SPECT imaging is limited to the solid collection angles defined by the collimators, generally, PET imaging has a higher sensitivity and spatial resolution than does SPECT. Therefore, the gold standard for spatial and time resolutions in nuclear imaging is defined for PET. For example, there is a gold standard for PET imaging for at rest myocardial perfusion with N-13-ammonia (NH3), at a dose of 740 MBq with attenuation correction.”
Conventional SPECT cameras generally employ an Anger camera, in which a single-pixel scintillation detector, such as NaI(Tl), LSO, GSO, CsI, CaF, or the like, is associated with a plurality of photomultipliers. Dedicated algorithms provide a two dimensional image of the scintillations in the single pixel scintillation detector. There are several disadvantages to this system, for example:    1. The dedicated algorithms associated with the single pixel cannot reach the accuracy of a two-dimensional image of a plurality of single pixel detectors;    2. The single-pixel detector is a rigid unit, which does not have the flexibility of motion of a plurality of small detectors, each with independent motion; and    3. A single hot spot may cause the single pixel detector of the Anger camera to saturate, whereas when a plurality of single pixel detectors is employed, saturation is localized to a few pixels and does not affect the whole image.
Other SPECT cameras which employ a plurality of single pixel detectors are also known.
U.S. Pat. No. 6,628,984, to Weinberg, issued on Sep. 30, 2003 and entitled, “Handheld camera with tomographic capability,” describes a tomographic imaging system, which includes a moveable detector or detectors capable of detecting gamma radiation; one or more position sensors for determining the position and angulation of the detector(s) in relation to a gamma ray emitting source; and a computational device for integrating the position and angulation of the detector(s) with information as to the energy and distribution of gamma rays detected by the detector and deriving a three dimensional representation of the source based on the integration. A method of imaging a radiation emitting lesion located in a volumetric region of interest also is disclosed.
U.S. Pat. No. 6,242,743, to DeVito, et al., issued on Jun. 5, 2001 and entitled, “Non-orbiting tomographic imaging system,” describes a tomographic imaging system which images ionizing radiation such as gamma rays or x rays and which: 1) can produce tomographic images without requiring an orbiting motion of the detector(s) or collimator(s) around the object of interest, 2) produces smaller tomographic systems with enhanced system mobility, and 3) is capable of observing the object of interest from sufficiently many directions to allow multiple time-sequenced tomographic images to be produced. The system consists of a plurality of detector modules which are distributed about or around the object of interest and which fully or partially encircle it. The detector modules are positioned close to the object of interest thereby improving spatial resolution and image quality. The plurality of detectors view a portion of the patient or object of interest simultaneously from a plurality of positions. These attributes are achieved by configuring small modular radiation detector with high-resolution collimators in a combination of application-specific acquisition geometries and non-orbital detector module motion sequences composed of tilting, swiveling and translating motions, and combinations of such motions. Various kinds of module geometry and module or collimator motion sequences are possible, and several combinations of such geometry and motion are shown. The geometric configurations may be fixed or variable during the acquisition or between acquisition intervals. Clinical applications of various embodiments of U.S. Pat. No. 6,242,743 include imaging of the human heart, breast, brain or limbs, or small animals. Methods of using the non-orbiting tomographic imaging system are also included.
U.S. Pat. No. 5,939,724, to Eisen, et al., issued on Aug. 17, 1999, and entitled, “Light weight-camera head and-camera assemblies containing it,” describes a lightweight gamma-camera head, assemblies, and kits that embody it. The gamma-camera head has a detector assembly which includes an array of room temperature, solid state spectroscopy grade detectors each associated with a collimator and preamplifier, which detectors and associated collimators and preamplifiers are arranged in parallel rows extending in a first direction and suitably spaced from each other in a second direction normal to the first direction, each of the parallel detector rows holding a plurality of detectors. The head may optionally have an electric motor for moving the detector in the second direction and optionally also in the first direction, either stepwise or continuously.
U.S. Pat. No. 6,525,320, to Juni, issued on Feb. 25, 2003, and entitled, single photon emission computed tomography system, describes a single photon emission computed tomography system, which produces multiple tomographic images of the type representing a three-dimensional distribution of a photon-emitting radioisotope. The system has a base including a patient support for supporting a patient such that a portion of the patient is located in a field of view. A longitudinal axis is defined through the field of view. A detector module is adjacent the field of view and includes a photon-responsive detector. The detector is an elongated strip with a central axis that is generally parallel to the longitudinal axis. The detector is operable to detect if a photon strikes the detector. The detector can also determine a position along the length of the strip where a photon is detected. A photon-blocking member is positioned between the field of view and the detector. The blocking member has an aperture slot for passage of photons aligned with the aperture slot. The slot is generally parallel to the longitudinal axis. A line of response is defined from the detector through the aperture. A displacement device moves either the detector module or the photon-blocking member relative to the other so that the aperture is displaced relative to the detector and the line of response is swept across at least a portion of the field of view.
U.S. Pat. No. 6,271,525, to Majewski, et al., issued on Aug. 7, 2001, and entitled, “Mini gamma camera, camera system and method of use,” describes a gamma camera, which comprises essentially and in order from the front outer or gamma ray impinging surface: 1) a collimator, 2) a scintillator layer, 3) a light guide, 4) an array of position sensitive, high resolution photomultiplier tubes, and 5) printed circuitry for receipt of the output of the photomultipliers. There is also described, a system wherein the output supplied by the high resolution, position sensitive photomultiplier tubes is communicated to: a) a digitizer and b) a computer where it is processed using advanced image processing techniques and a specific algorithm to calculate the center of gravity of any abnormality observed during imaging, and c) optional image display and telecommunications ports.
U.S. Pat. No. 6,271,524, to Wainer, et al., issued on Aug. 7, 2001 and entitled, “Gamma ray collimator,” describes a gamma ray collimator assembly comprising collimators of different gamma ray acceptance angles. For example, the acceptance angle of a first collimator may be between 0.2 and 5 degrees, and the acceptance angle of a second collimator may be between about 5 and 30 degrees.
U.S. Pat. No. 6,212,423, to Krakovitz, issued on Apr. 3, 2001 and entitled, “Diagnostic hybrid probes,” describes a hybrid nuclear and ultrasonic probe, comprising a cylindrical outer casing surrounding a nuclear probe, which comprises two scintillator plates intersecting perpendicularly, each of the scintillator plates having a plurality of parallel collimators; and an ultrasonic probe situated between said casing at the intersection of said scintillator plates.
List mode data acquisition is known in PET studies, and enables the determination of coincidence. It relates to recording every radiation event together with data pertinent to that event, which includes:    i. the time the radiation event impinged upon a detector pixel, with respect to a clock, with respect to a time bin, or with respect to another time definition, for example, a time interval between two clock signals; and    ii. the detector pixel location with respect to a coordinate system, at the time of the impinging.
The knowledge of time and location enables the determination of coincidence counts, namely photon counts that arrive substantially simultaneously, 180 degrees apart.
The time and location data may be stamped onto the radiation-event data packet, for example, as a header or as a footer, or otherwise associated with the radiation-event data packet, as known.
The time-stamped data available in PET studies may further be used for perfusion studies, where the timing of physiological processes of short durations, that is, durations shorter than about half the time span between heartbeats, is important. Perfusion studies usually involve a sequence of continuous acquisitions, each of which may represent data acquisition duration of about 10-30 seconds, although longer durations are sometimes employed. Data from each of the frames is independently reconstructed to form a set of images that can be visualized and used to estimate physiological parameters. This approach involves selection of the set of acquisition times, where one must choose between collecting longer scans with good counting statistics but poor temporal resolution, or shorter scans that are noisy but preserve temporal resolution.
US Patent Application 2003010539, to Tumer, et al., published on Jun. 5, 2003, and entitled, “X-ray and gamma ray detector readout system,” describes a readout electronics scheme, under development for high resolution, compact PET (positron emission tomography) imagers, using time tagging, based on LSO (lutetium ortho-oxysilicate, Lu.sub.2SiO.sub.5) scintillator and avalanche photodiode (APD) arrays.
There is some work relating to timing data in SPECT systems, employing Anger cameras.
U.S. Pat. No. 5,722,405, to Goldberg, issued on Mar. 3, 1998, and entitled, “Method and apparatus for acquisition and processing of event data in semi list mode,” describes a system for acquisition, processing and display of gated SPECT imaging data for use in diagnosing Coronary Artery Disease (CAD) in nuclear medicine, employing an Anger camera, and provides a physician with two parameters for evaluating CAD: information relating to the distribution of blood flow within the myocardium (perfusion) and information relating to myocardium wall motion (function). One aspect provides the physician with a display of functional images representing quantitative information relating to both perfusion and function with respect to selected regions of interest of the subject heart at end-diastole and end-systole segments of the cardiac cycle. The functional display consists of arcs of varied width depending on wall motion and color coded to illustrate degrees of myocardial perfusion for different pie shaped sections of a selected region of interest within a given short axis slice of reconstructed volumetric region data. Another aspect provides a series of display images allowing facilitated access, display, and comparison of the numerous image frames of the heart that may be collected during gated SPECT sessions. U.S. Pat. No. 5,722,405 also teaches the ability to define and recall parameter files representative of data acquisition and processing parameters and protocol for use in gated SPECT studies and includes a semi-list processing mode to increase efficiency of data acquisition within a camera computer system.
U.S. Pat. No. 7,026,623, to Oaknin, et al., issued on Apr. 11, 2006, and entitled, “Efficient single photon emission imaging,” describes-a method of diagnostic imaging in a shortened acquisition time for obtaining a reconstructed diagnostic image of a portion of a body of a human patient who has been administered with dosage of radiopharmaceutical substance radiating gamma rays, using SPECT and an Anger camera. The method comprises acquiring photons emitted from said portion of the body, by means of a detector capable of converting the photons into electric signals, wherein the total time of photon acquiring is substantially shorter than the clinically acceptable acquisition time; processing said electric signals by a position logic circuitry and thereby deriving data indicative of positions on said photon detector crystal, where the photons have impinged the detector; and reconstructing an image of a spatial distribution of the pharmaceutical substance within the portion of the body by iteratively processing said data. For example, the method includes effective acquisition time of less than 10 minutes, or less than 8 minutes, and acquiring photons in a list-mode procedure.
Current techniques record data with SPECT and electrocardiogram (ECG), and perform some gating to the data which is captured by the SPECT detectors, to incorporate the global and regional atrial and ventricular function and assessment of the relationship of perfusion to regional function.
Gated images are used to overcome distortions such as motion artifacts, which are caused due to motion of the heart during image acquisition. The images are needed as the physical model used for reconstruction assumes that the imaged objects are static. In gated imaging, photon-counting takes into account the portion of the heart contraction cycle within which a photon is measured. The Gating enables the reconstruction of an anatomical structure which is subject to periodic motion by enabling image acquisition only when the structure has reached the same configuration. Cardiac contraction is usually synchronized to the recorded electrocardiogram (ECG) signal that indicates the current heart pose. The period between a certain repetitive wave, such as R-wave, and a subsequent wave is divided into several time segments, called “frames”, which are usually spaced evenly. Each photon which is detected by the PET detectors during one of the frames is collected and associated with the related frame.
In gated imaging, each frame generates a single dataset. The collection of all the datasets belonging to all the frames are defined as a “dynamic” dataset.
The dynamic dataset is created by dividing the time span between one R-wave to the next R-wave into M frames that usually have an identical duration. Each detected photon is accumulated into a dataset of one of the M frames. Each dataset of the M datasets contains data relevant to a defined portion (“snapshot”) within the cardiac cycle.
Usually, during the image reconstruction process, each one of the gated datasets of the M frames is processed independently by a suitable reconstruction algorithm, see Leahy R et al., Computer tomography in: Handbook of Image and Video Processing, BovikA, Academic press, 2000, pp. 771-787; J. Kay. The EM algorithm in medical imaging, Stat. Meth. Med. Res., 6(1):55-75, January 1997; J. A. Fessler, Statistical image reconstruction methods for transmission tomography, Handbook of Medical Imaging, Volumetric region 2, pages 1-70. SPIE, Bellingham, 2000; R. M. Leahy et al., Statistical approaches in quantitative positron emission tomography, 10(2):14765, April 2000; M. Defrise, A short reader's guide to 3D tomographic reconstruction, Computerized Medical Imaging and Graphics, 25(2):1 13-6, March 2001; Vandenberghe, Y. D'Asseler, et al. Iterative reconstruction algorithms in nuclear medicine, Computerized Medical Imaging and Graphics, 25(2):105-11, March 2001; G. L. Zeng. Image reconstruction, a tutorial, Computerized Medical Imaging and Graphics, 25(2):97-103, March 2001; and R. M. Lewitt et al., Overview of methods for image reconstruction from projections in emission computed tomography, Proc. IEEE, 91(9):1588-611, October 2003, which are incorporated herein by reference in its entirety.
A common practice in gated SPECT reconstruction is to divide the gated dynamic dataset into M ‘non-gated’ data sets. Each one of the datasets includes data from a single frame i. The reconstruction of each volumetric region is performed independently using the relevant data set.
In particular, once the emission data is obtained, the data is processed to reconstruct the intensity distribution within the measured volumetric region. The reconstruction process is generally complex, due to the large quantity of data that must be processed in order to obtain an accurate reconstruction. The following prior art statistical model may be used to perform reconstruction.
We assume an intensity distribution, I, defined over an input overall volume U, where U denotes a set of basic elements, such as pixels in two dimensional overall volumes and voxels in three dimensional overall volumes, and I(u) is the intensity of a given basic element uεU. A detecting unit positioned on a radiation-emission-measuring-probe such as a PET detector or the like takes a series of measurements y=(yt)t=1T from different positions and orientations around the volumetric region U. The geometrical and physical properties of the detecting unit, together with its position and orientation in a given measurement t, determine the detection probability φt(u) of a photon emitted from location u in time t. Thus, the effective intensity of location u as viewed by the detecting unit during measurement t is φt(u)I(u).
The random count Xt(u) of photons that are emitted from location u and detected in measurement t is modeled by a Poisson process with mean φt(u)I(u). The total count of photons detected in measurement t is Yt=ΣuεU Xt(u), and the reconstruction problem is to reconstruct the intensities (I(u))uεU from the measurements (yt)t=1T.
The 2-D Radon transform is a mathematical relationship that may be used for reconstructing the emission intensities of volumetric region U when the set of measurements (yt)t=1T is unconstrained. The Radon transform is not statistical and does not take into account the Poissonian nature of the counts. In addition, it models the views as line projections. The Radon transform maps the spatial domain (x,y) to the Radon domain (p,φ). For a fixed projection angle, the Radon transform is simply a projection of the object. A technique known in the ART as filtered back-projection (FBP) uses a back-projection operator and the inverse of the Radon transform to reconstruct the intensity distribution in volumetric region U from measurements (yt)t=1T.
The basic, idealized problem solved by the FBP approach is to reconstruct an image from its Radon transform. The Radon transform, when properly defined, has a well-defined inverse. However, in order to invert the transform one needs measured data spanning 180°. In many medical imaging situations, the positioning of the detecting unit relative to the emitting object is constrained, so that complete measured data is not available. Reconstruction based on filtered back-projection is therefore of limited use for medical imaging. Maximum likelihood (ML) and Maximum A Posteriori (MAP) estimation methods, which address the statistical nature of the counts, have been found to provide better image reconstructions than FBP.
Limited-angle tomography is a reconstruction technique in the related art which reconstructs an image from projections acquired over a limited range of angular directions. The success of the reconstruction process depends upon the extent of the angular range acquired compared with the angular range of the missing projections. Any reconstruction from a limited range of projections potentially results in spatial distortions (artifacts) in the image. Limited angle techniques can be applied for both the Radon transform and the statistical models, but better results are generally achieved within the statistical framework. While it is known that the severity of the artifacts increases with the increasing angular range of the missing projections, limited-angle tomography does not provide information on which projections should be used in order to most effectively reconstruct the image.
ML estimation is a widely used method in the related art for reconstructing an image from a constrained set of measurements. A parameterization of the generative model described above is obtained by assigning an intensity I(u) to every voxel in U. The likelihood of the observed data y=(yt)t, given the set of parameters I={I(u):uεU} is:
                              L          ⁡                      (                          y              ❘              I                        )                          =                  ln          ⁢                                          ⁢                      P            ⁡                          (                              y                ❘                I                            )                                ⁢                                                                      =                                    ⁢                                                            ln                      ⁢                                                                                          ⁢                                                                        ∏                          t                                                ⁢                                                  P                          ⁡                                                      (                                                                                          y                                t                                                            ❘                              I                                                        )                                                                                                                =                                                                  ∑                        t                                            ⁢                                              ln                        ⁢                                                                                                  ⁢                                                  P                          ⁡                                                      (                                                                                                                            ∑                                  u                                                                ⁢                                                                                                      x                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                                              ❘                              I                                                        )                                                                                                                                                                                                                    =                                    ⁢                                                            ∑                      t                                        ⁢                                          ln                      ⁢                                                                                          ⁢                                              Poisson                        ⁡                                                  (                                                                                    y                              t                                                        ❘                                                                                          ∑                                u                                                            ⁢                                                                                                                                    ϕ                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                  ⁢                                                                  I                                  ⁡                                                                      (                                    u                                    )                                                                                                                                                                                )                                                                                                                                                                                          =                                    ⁢                                                            ∑                      t                                        ⁢                                          {                                                                                                                                                                  -                                                                                                      ∑                                    u                                                                    ⁢                                                                                                                                                    ϕ                                        t                                                                            ⁡                                                                              (                                        u                                        )                                                                                                              ⁢                                    I                                    ⁢                                                                          (                                      u                                      )                                                                                                                                                                  +                                                                                                                                                                                                                                                            y                                  t                                                                ⁢                                ln                                ⁢                                                                                                                                  ⁢                                                                                                      ∑                                    u                                                                    ⁢                                                                                                                                                    ϕ                                        t                                                                            ⁡                                                                              (                                        u                                        )                                                                                                              ⁢                                                                          I                                      ⁡                                                                              (                                        u                                        )                                                                                                                                                                                                        -                                                              ln                                ⁡                                                                  (                                                                                                            y                                      t                                                                        !                                                                    )                                                                                                                                                                                        }                                                                                                                              (        1        )            
Note that the lower and upper bound of an indexing variable, such as voxels u and time index t, are omitted in the following description, when they are clear from the context.
There is currently no analytic way to solve Eqn. 1 for the maximum of the likelihood function. However, optimization methods that find local maxima of the likelihood are known. One such method is the Expectation-Maximization (EM) process.
Since the data generated by the model is only partially observable by our measurements, a basic ingredient of the EM formalism is to define a set of random variables that completely define the data generated by the model. In the current case, since Yt=ΣuXt(u), the set of variables {Xu(t):uεU; t=1, . . . , T} is such a set; the generated data is x=(xt)t, where xt=(xt(u))u, and the observed data y is completely determined by x. The main tool in the EM formalism is the complete data likelihood:
                                                                                             ln                  ⁢                                                                          ⁢                                      P                    ⁡                                          (                                              x                        ❘                        I                                            )                                                                      =                                ⁢                                                      ln                    ⁢                                                                  ∏                        t                                            ⁢                                              P                        ⁡                                                  (                                                                                    x                              t                                                        ❘                            I                                                    )                                                                                                      =                                                            ∑                      t                                        ⁢                                          ln                      ⁢                                                                        ∏                          u                                                ⁢                                                  Poisson                          ⁡                                                      (                                                                                                                            x                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ❘                                                                                                                                    ϕ                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                  ⁢                                                                  I                                  ⁡                                                                      (                                    u                                    )                                                                                                                                                        )                                                                                                                                                                                                                      =                                ⁢                                                      ∑                    t                                    ⁢                                                            ∑                      u                                        ⁢                                          {                                                                                                    -                                                                                          ϕ                                t                                                            ⁡                                                              (                                u                                )                                                                                                              ⁢                                                      I                            ⁡                                                          (                              u                              )                                                                                                      +                                                                                                            x                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                      ln                            ⁡                                                          (                                                                                                                                    ϕ                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                  ⁢                                                                  I                                  ⁡                                                                      (                                    u                                    )                                                                                                                              )                                                                                                      +                                                  ln                          ⁡                                                      (                                                                                                                            x                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            !                                                        )                                                                                              }                                                                                                                                (          2          )                    
Since the likelihood depends on the complete data, which is only partially observable, we take its expectation with respect to the overall volume of the unobserved data, given the current set of hypothesized parameters (i.e. the current estimator). The result is a function Q(I|I′) which assigns likelihood to sets I of model parameters, given the current set I′, and given the observed data y:
                                                                                             Q                  ⁡                                      (                                          I                      ❘                                              I                        ′                                                              )                                                  =                                ⁢                                  E                  ⁡                                      [                                                                                            ln                          ⁢                                                                                                          ⁢                                                      P                            ⁡                                                          (                                                              x                                ❘                                I                                                            )                                                                                                      ❘                        y                                            ;                                              I                        ′                                                              ]                                                                                                                          =                                ⁢                                                      ∑                    t                                    ⁢                                                            ∑                      u                                        ⁢                                          {                                                                                                                                                                                                    -                                                                                                            ϕ                                      t                                                                        ⁡                                                                          (                                      u                                      )                                                                                                                                      ⁢                                                                  I                                  ⁡                                                                      (                                    u                                    )                                                                                                                              +                                                                                                                                                                                                                                                            E                                  ⁡                                                                      [                                                                                                                                                                                                        x                                            t                                                                                    ⁡                                                                                      (                                            u                                            )                                                                                                                          ❘                                                                                  y                                          t                                                                                                                    ;                                                                              I                                        ′                                                                                                              ]                                                                                                  ⁢                                                                  ln                                  ⁡                                                                      (                                                                                                                                                            ϕ                                          t                                                                                ⁡                                                                                  (                                          u                                          )                                                                                                                    ⁢                                                                              I                                        ⁡                                                                                  (                                          u                                          )                                                                                                                                                      )                                                                                                                              +                              C                                                                                                                          }                                                                                                                                (          3          )                    
where C is a term which is independent of the intensities I. The function Q(I|I′) is maximized by the following new estimates:
                                          I            ⁡                          (              u              )                                =                                    1                                                ∑                  t                                ⁢                                                      ϕ                    t                                    ⁡                                      (                    u                    )                                                                        ⁢                                          ∑                t                            ⁢                              E                ⁡                                  [                                                                                                              x                          t                                                ⁡                                                  (                          u                          )                                                                    ❘                                              y                        t                                                              ;                                          I                      ′                                                        ]                                                                    ;                  ∀                      u            ∈                          U              .                                                          (        4        )            
The expectation in Equation 4 is obtained as follows:
                                                                                                                 P                                                                  X                        t                                            ⁡                                              (                        u                        )                                                                              ⁡                                      (                                                                                                                        x                            t                                                    ⁡                                                      (                            u                            )                                                                          ❘                                                  y                          t                                                                    ;                                              I                        ′                                                              )                                                  =                                ⁢                                                                                                    P                                                  Y                          t                                                                    ⁡                                              (                                                                                                            y                              t                                                        ❘                                                                                          x                                t                                                            ⁡                                                              (                                u                                )                                                                                                              ;                                                      I                            ′                                                                          )                                                              ⁢                                                                  P                                                                              X                            t                                                    ⁡                                                      (                            u                            )                                                                                              ⁡                                              (                                                                                                            x                              t                                                        ⁡                                                          (                              u                              )                                                                                ❘                                                      I                            ′                                                                          )                                                                                                                        P                                              Y                        t                                                              ⁡                                          (                                                                        y                          t                                                ❘                                                  I                          ′                                                                    )                                                                                                                                              =                                ⁢                                                                                                                              Poisson                          (                                                                                                                    y                                t                                                            -                                                                                                x                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                                                      ❘                                                                                          ∑                                                                  v                                  ≠                                  u                                                                                            ⁢                                                                                                                                    ϕ                                    t                                                                    ⁡                                                                      (                                    v                                    )                                                                                                  ⁢                                                                                                      I                                    ′                                                                    ⁡                                                                      (                                    v                                    )                                                                                                                                                                                )                                                                                                                                                              Poisson                          ⁢                                                                                                          ⁢                                                      (                                                                                                                            x                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ❘                                                                                                                                    ϕ                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                  ⁢                                                                                                      I                                    ′                                                                    ⁡                                                                      (                                    u                                    )                                                                                                                                                        )                                                                                                                                                    Poisson                    ⁢                                                                                  ⁢                                          (                                                                        y                          t                                                ❘                                                                              ∑                            v                                                    ⁢                                                                                                                    ϕ                                t                                                            ⁡                                                              (                                v                                )                                                                                      ⁢                                                          I                              ⁡                                                              (                                v                                )                                                                                                                                                        )                                                                                                                                              =                                ⁢                                  Binomial                  ⁡                                      (                                                                                                                        x                            t                                                    ⁡                                                      (                            u                            )                                                                          ❘                                                                                                                                            ϕ                                t                                                            ⁡                                                              (                                u                                )                                                                                      ⁢                                                                                          I                                ′                                                            ⁡                                                              (                                u                                )                                                                                                                                                                        ∑                              v                                                        ⁢                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  v                                  )                                                                                            ⁢                                                                                                I                                  ′                                                                ⁡                                                                  (                                  v                                  )                                                                                                                                                                                        ;                                              y                        t                                                              )                                                                                                            (          5          )                    
It follows that
            E      ⁡              [                                            xt              ⁡                              (                u                )                                      ❘            yt                    ;                      I            ′                          ]              =                  y        t            ⁢                                                  ϕ              t                        ⁡                          (              u              )                                ⁢                                    I              ′                        ⁡                          (              u              )                                                            ∑            v                    ⁢                                                    ϕ                t                            ⁡                              (                v                )                                      ⁢                                          I                ′                            ⁡                              (                v                )                                                          ,and hence the EM iteration is:
                              I          ⁡                      (            u            )                          =                              1                                          ∑                t                            ⁢                                                ϕ                  t                                ⁡                                  (                  u                  )                                                              ⁢                                    ∑              t                        ⁢                                          y                t                            ⁢                                                                                          ϕ                      t                                        ⁡                                          (                      u                      )                                                        ⁢                                                            I                      ′                                        ⁡                                          (                      u                      )                                                                                                            ∑                    v                                    ⁢                                                                                    ϕ                        t                                            ⁡                                              (                        v                        )                                                              ⁢                                                                  I                        ′                                            ⁡                                              (                        v                        )                                                                                                                                                    (        6        )            
It is provable that each EM iteration improves the likelihood. Thus, given a random starting estimator, the EM algorithm iterates the above improvement step until it converges to a local maximum of the likelihood. Several random starts increase the chance of finding a globally good estimator.
It is usually desired to maximize the expected posterior probability (given a proper prior) rather than the expected likelihood. In that case we assume a prior probability on the intensities P(I)=ΠuP(I(u)). A proper conjugate prior for the Poisson distribution is the Gamma distribution:
                              P          ⁡                      (                          I              ⁡                              (                u                )                                      )                          =                              Gamma            ⁢                                                  ⁢                          (                                                                    I                    ⁡                                          (                      u                      )                                                        ❘                                      α                    u                                                  ;                                  β                  u                                            )                                =                                                    β                u                                                      α                    u                                    +                  1                                                            Γ                ⁡                                  (                                                            α                      u                                        +                    1                                    )                                                      ⁢                                          I                ⁡                                  (                  u                  )                                                            α                u                                      ⁢                          ⅇ                                                -                                      β                    u                                                  ⁢                                  I                  ⁡                                      (                    u                    )                                                                                                          (        7        )            
Now the maximization is done on Q(I|I′)=E[lnP(x|I)p(I)|y; I′]. Plugging the Gamma prior into Q, and solving for I(u), we get the following EM iteration for the maximum posterior estimation:
                              I          ⁡                      (            u            )                          =                                            α              u                        +                                          ∑                t                            ⁢                              E                ⁡                                  [                                                                                                              x                          t                                                ⁡                                                  (                          u                          )                                                                    ❘                                              y                        t                                                              ;                                          I                      ′                                                        ]                                                                                        β              u                        +                          ∑                                                ϕ                  t                                ⁡                                  (                  u                  )                                                                                        (        8        )                                =                              1                                          β                u                            +                                                ∑                  t                                ⁢                                                      ϕ                    t                                    ⁡                                      (                    u                    )                                                                                ⁡                      [                                          α                u                            +                                                ∑                  t                                ⁢                                                      y                    t                                    ⁢                                                                                                              ϕ                          t                                                ⁡                                                  (                          u                          )                                                                    ⁢                                                                        I                          ′                                                ⁡                                                  (                          u                          )                                                                                                                                    ∑                        v                                            ⁢                                                                                                    ϕ                            t                                                    ⁡                                                      (                            u                            )                                                                          ⁢                                                                              I                            ′                                                    ⁡                                                      (                            v                            )                                                                                                                                                                    ]                                              (        9        )            
The EM update step can be formulated in matrix notation as follows. Let Φ be the matrix of the projections [φt(u)]t, u, and let I, I′, y, α and β be represented as column vectors. Equation 8 can be written in vector and matrix notations as:
                    I        =                              α            +                                          I                ′                            ·                              (                                                      Φ                    T                                    ⁢                                      y                                          Φ                      ⁢                                                                                          ⁢                                              I                        ′                                                                                            )                                                          β            +                                          Φ                T                            ⁢              1                                                          (        10        )            
where the explicit multiplication and division denote element-wise operations, and where 1 is a vector (of the appropriate length) consisting solely of 1's.
Limited computational resources (i.e., when the entire projection matrix Φ cannot be kept in memory) may require breaking the update computation according to a partition of Φ into a set of sub-matrices (Φi). In that case the intensities can be updated gradually (using only one sub-matrix at each step) according to the following computation:
                    I        =                              α            +                                          I                ′                            ·                                                ∑                  i                                ⁢                                                      Φ                    i                    T                                    ⁢                                                            y                      t                                                                                      Φ                        i                                            ⁢                                              I                        ′                                                                                                                                      β            +                                          ∑                i                            ⁢                                                Φ                  i                  T                                ⁢                1                                                                        (        11        )            
where yi denotes the vector of observations that are obtained using the views of Φi.
In order to achieve a reconstructed image which is adequate for medical diagnostic and treatment purposes, a high-resolution image of the tested object must be obtained. When high-resolution detecting units are used, their efficiency is relatively low, and the detecting units must remain at each position for a relatively long time in order to achieve a high probability of detection. Since during medical testing, measurements are generally performed at many locations as the detecting unit is moved relative to the observed organ, the testing procedure generally requires a long time and is physically and emotionally difficult for the patient. Additionally, reconstruction is based upon a large quantity of data, and is a lengthy and computationally complex process.
Reference is now made to FIG. 13, which is a schematic flowchart that illustrates steps of a typical prior art gated image reconstruction method. In FIG. 13, i denotes a frame counter and M denotes the number of frames. Usually, after all the M datasets are fetched, as shown at 2, and i is set with 1, as shown at 4, the dataset that corresponds with frame i is loaded into the processing unit, as shown at 6. During the following step, as shown at 8, the processing unit is used to perform a frame reconstruction according to the dataset that corresponds with frame i using any suitable reconstruction algorithm known in the art, such as the aforementioned EM algorithm, ordered subset expectation maximization (OSEM), and algebraic reconstruction techniques (ART) FBP.
As shown at 12, after the reconstruction is completed, the frame counter i is incremented by 1. If i is larger than M and there are no more frame the process ends. If i is not larger than M, the next dataset that corresponds with the subsequent frame is loading for reconstruction. In such a manner, the generation of a static imaging of the heart in a specific configuration becomes possible.
However, a known problem of such a method is the high computational load that is needed for the execution thereof. Even when using an ordered set method, such as the aforementioned Hudson et al. method, which is capable of reducing the computational load while giving similar results, the computational power needed to obtain a good image reconstruction is still quite high. Such a high computational load results in a longer reconstruction time that reduces the throughput of the processing unit and requires a more expensive processing hardware.
There is thus a widely recognized need for, and it would be highly advantageous to have, a method and an apparatus for image reconstruction in nuclear medicine imaging devoid of the above limitations.